{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mRecompiling stale cache file /home/ehsan/.julia/lib/v0.6/Conda.ji for module Conda.\n",
      "\u001b[39m\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mModule Compat with uuid 2711308643704 is missing from the cache.\n",
      "This may mean module Compat does not support precompilation but is imported by a module that does.\u001b[39m\n",
      "\u001b[1m\u001b[91mERROR: \u001b[39m\u001b[22mLoadError: \u001b[91mDeclaring __precompile__(false) is not allowed in files that are being precompiled.\u001b[39m\n",
      "Stacktrace:\n",
      " [1] \u001b[1m_require\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Symbol\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:448\u001b[22m\u001b[22m\n",
      " [2] \u001b[1mrequire\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Symbol\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:398\u001b[22m\u001b[22m\n",
      " [3] \u001b[1minclude_from_node1\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:569\u001b[22m\u001b[22m\n",
      " [4] \u001b[1minclude\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./sysimg.jl:14\u001b[22m\u001b[22m\n",
      " [5] \u001b[1manonymous\u001b[22m\u001b[22m at \u001b[1m./<missing>:2\u001b[22m\u001b[22m\n",
      "while loading /home/ehsan/.julia/v0.6/Conda/src/Conda.jl, in expression starting on line 34\n"
     ]
    },
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mFailed to precompile Conda to /home/ehsan/.julia/lib/v0.6/Conda.ji.\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mFailed to precompile Conda to /home/ehsan/.julia/lib/v0.6/Conda.ji.\u001b[39m",
      "",
      "Stacktrace:",
      " [1] \u001b[1mcompilecache\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:703\u001b[22m\u001b[22m",
      " [2] \u001b[1m_require\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Symbol\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:456\u001b[22m\u001b[22m",
      " [3] \u001b[1mrequire\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Symbol\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:398\u001b[22m\u001b[22m",
      " [4] \u001b[1m_include_from_serialized\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:157\u001b[22m\u001b[22m",
      " [5] \u001b[1m_require_from_serialized\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Int64, ::Symbol, ::String, ::Bool\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:200\u001b[22m\u001b[22m",
      " [6] \u001b[1m_require_search_from_serialized\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Int64, ::Symbol, ::String, ::Bool\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:236\u001b[22m\u001b[22m",
      " [7] \u001b[1m_require\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Symbol\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:434\u001b[22m\u001b[22m",
      " [8] \u001b[1mrequire\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Symbol\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:398\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "using PyPlot, Polynomials, CoolProp, Roots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thermodynamics\n",
    "This section finds the correlations for the density, viscosity, and partition coefficient (in water/oil) of DME.  \n",
    "### Note: \n",
    "All the data in the Shell paper is reported as mole fraction. For convenience, I convert them to mass fraction, assuming the molecular weight of crude oil is equal to the molecular weight of C10H22."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A detour: viscosity measurements\n",
    "Is it possible to measure DME-water viscosity in a tube?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "flow rate is 9.42477796076938 ml/min\n",
      "pressure drop is 63.99999999999999 mbar\n"
     ]
    }
   ],
   "source": [
    "Re = 200 # Reynolds number\n",
    "ρ  = 1000 # [kg/m^3] density\n",
    "L = 1.0   # [m] pipe length\n",
    "D = 0.001 # [m] pipe diameter\n",
    "μ = 0.001 # [Pa.s]\n",
    "v = Re*μ/(ρ*D) # velocity\n",
    "q = v*π*D^2/4  # [m^3/s] flow rate\n",
    "println(\"flow rate is $(q*1e6*60) ml/min\")\n",
    "dp = 8*μ*L*q/(π*(D/2)^4) # [Pa] pressure drop\n",
    "println(\"pressure drop is $(dp/100) mbar\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.04607"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# some data\n",
    "MW_water = 0.018 # [kg/mol]\n",
    "MW_oil   = 0.142 # [kg/mol]   \n",
    "MW_DME   = 0.04607 # [kg/mol]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Density of DME-water and DME-oil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: polyfit not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: polyfit not defined\u001b[39m",
      ""
     ]
    }
   ],
   "source": [
    "# Fig 4, SPE-177919, T = 70 degC, p = 2000 psi\n",
    "x_DME_oil = [0.0, 0.04, 0.10, 0.39, 0.62]         # [mol frac]\n",
    "w_DME_oil = MW_DME*x_DME_oil./(MW_DME*x_DME_oil+MW_oil*(1-x_DME_oil))\n",
    "ρ_oil     = [0.748, 0.746, 0.744, 0.722, 0.697]*1000   # [g/cm^3]*1000 = [kg/m^3]\n",
    "\n",
    "rho_fit = polyfit(w_DME_oil, ρ_oil, 2)\n",
    "\n",
    "plot(w_DME_oil, ρ_oil, \"o\", w_DME_oil, polyval(rho_fit, w_DME_oil))\n",
    "xlabel(\"DME mass fraction in oil [-]\")\n",
    "ylabel(\"Oil density [g/cm^3]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## K-value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: polyfit not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: polyfit not defined\u001b[39m",
      ""
     ]
    }
   ],
   "source": [
    "# experimental data from SPE-179771\n",
    "# data for 1.4% brine and 71 degC and 2000 psia\n",
    "x_DME_oil = [0.15, 0.30, 0.5, 0.6]   # [mol frac] [0.5, 13] is a madeup point\n",
    "w_DME_oil = MW_DME*x_DME_oil./(MW_DME*x_DME_oil+MW_oil*(1-x_DME_oil)) # [mass frac]\n",
    "K_value   = [17.0, 15.0, 13.0, 12.5]  # [-]\n",
    "x_DME_water = x_DME_oil./K_value # [mole frac]\n",
    "w_DME_water = MW_DME*x_DME_water./(MW_DME*x_DME_water+MW_water*(1-x_DME_oil))\n",
    "K_value_mass = w_DME_oil./w_DME_water\n",
    "K_fit = polyfit(w_DME_water, K_value_mass, 2)\n",
    "w_DME_plot = collect(linspace(minimum(w_DME_water), maximum(w_DME_water), 50))\n",
    "plot(w_DME_water, K_value_mass, \"o\", w_DME_plot, polyval(K_fit,w_DME_plot))\n",
    "xlabel(\"DME mass fraction in water [-]\")\n",
    "ylabel(\"K-value [-]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# viscosity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: polyfit not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: polyfit not defined\u001b[39m",
      ""
     ]
    }
   ],
   "source": [
    "# Fig 4, SPE-177919, T = 70 degC, p = 2000 psi\n",
    "x_DME_oil = [0.0, 0.04, 0.10, 0.39, 0.62] # [mol frac]\n",
    "w_DME_oil = MW_DME*x_DME_oil./(MW_DME*x_DME_oil+MW_oil*(1-x_DME_oil))\n",
    "μ_oil = [1.04, 0.96, 0.84, 0.52, 0.34]/1000.0    # [cP]/1000 = [Pa.s]\n",
    "mu_oil_fit = polyfit(w_DME_oil, μ_oil, 2)\n",
    "w_DME_plot = collect(linspace(minimum(w_DME_oil), maximum(w_DME_oil), 50))\n",
    "plot(w_DME_oil, μ_oil, \"o\", w_DME_plot, polyval(mu_oil_fit, w_DME_plot))\n",
    "xlabel(\"DME mass fraction in oil [-]\")\n",
    "ylabel(\"Oil viscosity [Pa.s]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the viscosity of water with DME, I use the package coolprop to calculate the water viscosity at 70 degC and 2000 psi and assume a linear increase of viscosity when water is mixed with DME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: PropsSI not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: PropsSI not defined\u001b[39m",
      ""
     ]
    }
   ],
   "source": [
    "T0 = 70 + 273.15 # [K] \n",
    "p0 = 2000/14.7*1e5   # [Pa]\n",
    "mu_water = PropsSI(\"viscosity\", \"T\", T0, \"P\", p0, \"water\") # [Pa.s]\n",
    "mu_water_DME = [mu_water, 2*mu_water]\n",
    "w_DME_water = [0.0, maximum(w_DME_oil)/maximum(K_value_mass)]\n",
    "mu_water_fit = polyfit(w_DME_water, mu_water_DME, 1)\n",
    "plot(w_DME_water, mu_water_DME, \"o\", w_DME_water, mu_water_fit(w_DME_water))\n",
    "xlabel(\"DME mass fraction in water [-]\")\n",
    "ylabel(\"Water-DME viscosity [Pa.s]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The PDE's\n",
    "All the PDE's are explained in the documents. In this section, I use the JFVM package to discretize them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPrecompiling module Conda.\n",
      "\u001b[39m\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mModule Compat with uuid 2711308643704 is missing from the cache.\n",
      "This may mean module Compat does not support precompilation but is imported by a module that does.\u001b[39m\n",
      "\u001b[1m\u001b[91mERROR: \u001b[39m\u001b[22mLoadError: \u001b[91mDeclaring __precompile__(false) is not allowed in files that are being precompiled.\u001b[39m\n",
      "Stacktrace:\n",
      " [1] \u001b[1m_require\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Symbol\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:448\u001b[22m\u001b[22m\n",
      " [2] \u001b[1mrequire\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Symbol\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:398\u001b[22m\u001b[22m\n",
      " [3] \u001b[1minclude_from_node1\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:569\u001b[22m\u001b[22m\n",
      " [4] \u001b[1minclude\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./sysimg.jl:14\u001b[22m\u001b[22m\n",
      " [5] \u001b[1manonymous\u001b[22m\u001b[22m at \u001b[1m./<missing>:2\u001b[22m\u001b[22m\n",
      "while loading /home/ehsan/.julia/v0.6/Conda/src/Conda.jl, in expression starting on line 34\n"
     ]
    },
    {
     "ename": "LoadError",
     "evalue": "LoadError: \u001b[91mFailed to precompile Conda to /home/ehsan/.julia/lib/v0.6/Conda.ji.\u001b[39m\nwhile loading /home/ehsan/.julia/v0.6/JFVM/src/JFVM.jl, in expression starting on line 3",
     "output_type": "error",
     "traceback": [
      "LoadError: \u001b[91mFailed to precompile Conda to /home/ehsan/.julia/lib/v0.6/Conda.ji.\u001b[39m\nwhile loading /home/ehsan/.julia/v0.6/JFVM/src/JFVM.jl, in expression starting on line 3",
      "",
      "Stacktrace:",
      " [1] \u001b[1mcompilecache\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:703\u001b[22m\u001b[22m",
      " [2] \u001b[1m_require\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Symbol\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:490\u001b[22m\u001b[22m",
      " [3] \u001b[1mrequire\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Symbol\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:398\u001b[22m\u001b[22m",
      " [4] \u001b[1m_include_from_serialized\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:157\u001b[22m\u001b[22m",
      " [5] \u001b[1m_require_from_serialized\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Int64, ::Symbol, ::String, ::Bool\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:200\u001b[22m\u001b[22m",
      " [6] \u001b[1m_require_search_from_serialized\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Int64, ::Symbol, ::String, ::Bool\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:236\u001b[22m\u001b[22m",
      " [7] \u001b[1m_require\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Symbol\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:434\u001b[22m\u001b[22m",
      " [8] \u001b[1mrequire\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Symbol\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:398\u001b[22m\u001b[22m",
      " [9] \u001b[1minclude_from_node1\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:569\u001b[22m\u001b[22m",
      " [10] \u001b[1meval\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Module, ::Any\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./boot.jl:235\u001b[22m\u001b[22m",
      " [11] \u001b[1m_require\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Symbol\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:483\u001b[22m\u001b[22m",
      " [12] \u001b[1mrequire\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Symbol\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:398\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "using JFVM\n",
    "include(\"../rel_perms_real.jl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: PropsSI not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: PropsSI not defined\u001b[39m",
      ""
     ]
    }
   ],
   "source": [
    "# T = 70 + 273.15 # [K] \n",
    "# p = 2000/14.7*1e5   # [Pa]\n",
    "rho_water = PropsSI(\"D\", \"T\", T0, \"P\", p0, \"water\")\n",
    "# Corey rel-perm parameters\n",
    "krw0 = 0.2\n",
    "kro0 = 0.8\n",
    "n_o  = 2.0\n",
    "n_w  = 2.0\n",
    "swc  = 0.08\n",
    "sor  = 0.3\n",
    "sorting_factor = 2.4\n",
    "pce = 100 # [Pa]\n",
    "pc0 = 1e5 # [Pa]\n",
    "contact_angle = deg2rad(20) # [radian]\n",
    "\n",
    "perm_val  = 0.01e-12 # [m^2] permeability\n",
    "poros_val = 0.40     # [-] porosity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(::#33) (generic function with 1 method)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "KRW  = sw -> krw(sw, krw0, sor, swc, n_w)\n",
    "dKRW = sw -> dkrwdsw(sw, krw0, sor, swc, n_w)\n",
    "KRO  = sw -> kro(sw, kro0, sor, swc, n_o)\n",
    "dKRO = sw -> dkrodsw(sw, kro0, sor, swc, n_o)\n",
    "# PC   = sw -> pc_imb(sw, pce, swc, sor, teta=contact_angle, \n",
    "#                 labda=sorting_factor, b = 0.0, pc01=pc0, pc02=pc0)\n",
    "# dPC  = sw -> dpc_imb(sw, pce, swc, sor, teta=contact_angle, \n",
    "#                 labda1=sorting_factor, labda2=sorting_factor,\n",
    "#                 b = 1.0, pc01=pc0, pc02=pc0)\n",
    "# PC2  = sw -> pc_imb2(sw, pce, swc, sor, teta=contact_angle, \n",
    "#                 labda=sorting_factor, b = 0.0, pc_star1=pc0, \n",
    "#                 pc_star2=pc0)\n",
    "# dPC2 = sw -> dpc_imb2(sw, pce, swc, sor, teta=contact_angle, \n",
    "#                 labda1=sorting_factor, labda2=sorting_factor, \n",
    "#                 b = 1.0, pc_star1=pc0, pc_star2=pc0)\n",
    "PC  = sw -> pc_imb3(sw, pce, swc, sor, teta=contact_angle, \n",
    "                labda=sorting_factor, b = 0.0, pc_star1=pc0, \n",
    "                pc_star2=pc0)\n",
    "dPC = sw -> dpc_imb3(sw, pce, swc, sor, teta=contact_angle, \n",
    "                labda1=sorting_factor, labda2=sorting_factor, \n",
    "                b = 1.0, pc_star1=pc0, pc_star2=pc0)\n",
    "d2PC = sw -> d2pc_imb3(sw, pce, swc, sor, teta=contact_angle, \n",
    "                labda1=sorting_factor, labda2=sorting_factor, \n",
    "                b = 1.0, pc_star1=pc0, pc_star2=pc0)\n",
    "# PCdrain2 = sw -> pc_drain2(sw, pce, swc, labda=sorting_factor, pc_star=pc0)\n",
    "# PCdrain3 = sw -> pc_drain3(sw, pce, swc, labda=sorting_factor, pc_star=pc0)\n",
    "ρ_oil = w_DME_oil -> polyval(rho_fit, w_DME_oil)\n",
    "dρ_oil = w_DME_oil -> polyval(polyder(rho_fit), w_DME_oil)\n",
    "μ_water = w_DME_water -> polyval(mu_water_fit, w_DME_water)\n",
    "dμ_water = w_DME_water -> polyval(polyder(mu_water_fit), w_DME_water)\n",
    "μ_oil = w_DME_oil -> polyval(mu_oil_fit, w_DME_oil)\n",
    "dμ_oil = w_DME_oil -> polyval(polyder(mu_oil_fit), w_DME_oil)\n",
    "Kval = w_DME_water -> polyval(K_fit, w_DME_water)\n",
    "dKval = w_DME_water -> polyval(polyder(K_fit), w_DME_water)\n",
    "w_oil = w_DME_water -> polyval(Poly([0.0,1.0])*K_fit, w_DME_water)\n",
    "dw_oil = w_DME_water -> polyval(polyder(Poly([0.0,1.0])*K_fit), w_DME_water)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: Poly not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: Poly not defined\u001b[39m",
      ""
     ]
    }
   ],
   "source": [
    "Poly([0.0,1.0])*Poly([1.0,1.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: plot not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: plot not defined\u001b[39m",
      ""
     ]
    }
   ],
   "source": [
    "sw_plot = collect(linspace(0,1,1000))\n",
    "plot(sw_plot, KRW.(sw_plot), sw_plot, KRO.(sw_plot))\n",
    "axis([0,1,0,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: plot not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: plot not defined\u001b[39m",
      ""
     ]
    }
   ],
   "source": [
    "plot(sw_plot, PC.(sw_plot))\n",
    "plot([swc, swc], [-pc0, pc0], [1-sor, 1-sor], [-pc0, pc0], \n",
    "    [0,1], [0,0], \"--\")\n",
    "xlabel(\"Water saturation (Sw) [-]\")\n",
    "ylabel(\"Capillary pressure (Pc) [Pa]\")\n",
    "sw_imb_end = fzero(PC, [swc, 1-sor])\n",
    "plot(sw_imb_end, 0, \"o\")\n",
    "# axis([0, 1, -pc0, pc0])\n",
    "yscale(\"symlog\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: createMesh1D not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: createMesh1D not defined\u001b[39m",
      ""
     ]
    }
   ],
   "source": [
    "Lx   = 1.0 # [m]\n",
    "Nx  = 20  # number of grids in the x direction\n",
    "m   = createMesh1D(Nx, Lx)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boundary conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: createBC not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: createBC not defined\u001b[39m",
      ""
     ]
    }
   ],
   "source": [
    "u_inj = 1.0/(3600*24) # 1 m/day to m/s injection velocity\n",
    "c_inj = 0.2 # DME mass fraction in the injected water\n",
    "\n",
    "BCp = createBC(m) # pressure boundary condition\n",
    "BCs = createBC(m) # saturation boundary condition\n",
    "BCc = createBC(m) # concentration (DME in water) boundary condition\n",
    "\n",
    "BCp.right.a[:] = 0.0\n",
    "BCp.right.b[:] = 1.0\n",
    "BCp.right.c[:] = p0\n",
    "\n",
    "BCp.left.a[:]  = perm_val/μ_water(c_inj)\n",
    "BCp.left.b[:]  = 0.0\n",
    "BCp.left.c[:]  = -u_inj\n",
    "\n",
    "BCs.left.a[:]  = 0.0\n",
    "BCs.left.b[:]  = 1.0\n",
    "BCs.left.c[:]  = 1.0\n",
    "\n",
    "BCs.right.a[:]  = 0.0\n",
    "BCs.right.b[:]  = 1.0\n",
    "BCs.right.c[:]  = sw_imb_end\n",
    "\n",
    "BCc.left.a[:]  = 0.0\n",
    "BCc.left.a[:]  = 1.0\n",
    "BCc.left.a[:]  = c_inj\n",
    "\n",
    "# discretize\n",
    "M_bc_p, RHS_bc_p = boundaryConditionTerm(BCp)\n",
    "M_bc_s, RHS_bc_s = boundaryConditionTerm(BCs)\n",
    "M_bc_c, RHS_bc_c = boundaryConditionTerm(BCc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: swc not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: swc not defined\u001b[39m",
      ""
     ]
    }
   ],
   "source": [
    "c0  = 0.0 # [mass frac] initial DME concentration in the water phase\n",
    "sw0 = swc+0.05 # [vol frac] initial water saturation\n",
    "\n",
    "p_init  = createCellVariable(m, p0)\n",
    "c_init  = createCellVariable(m, c0)\n",
    "sw_init = createCellVariable(m, sw0)\n",
    "\n",
    "# new values of each variable\n",
    "p_val  = createCellVariable(m, p0)\n",
    "c_val  = createCellVariable(m, c0)\n",
    "sw_val = createCellVariable(m, sw0)\n",
    "\n",
    "Δsw_init = sw_init - sw_val # we solve for this variable\n",
    "Δc_init  = c_init - c_val   # we solve for this variable\n",
    "\n",
    "# Other cell variables\n",
    "k = createCellVariable(m, perm_val)\n",
    "ϕ = createCellVariable(m, poros_val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solver parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: poros_val not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: poros_val not defined\u001b[39m",
      ""
     ]
    }
   ],
   "source": [
    "n_pv    = 2.0 # number of injected pore volumes\n",
    "t_final = n_pv*Lx/(u_inj/poros_val) # [s] final time\n",
    "dt      = t_final/n_pv/Nx # [s] time step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All the equation that I need for the linearization\n",
    "ρ_oil = w_DME_oil -> polyval(rho_fit, w_DME_oil)\n",
    "dρ_oil = w_DME_oil -> polyval(polyder(rho_fit), w_DME_oil)\n",
    "μ_water = w_DME_water -> polyval(mu_water_fit, w_DME_water)\n",
    "dμ_water = w_DME_water -> polyval(polyder(mu_water_fit), w_DME_water)\n",
    "μ_oil = w_DME_oil -> polyval(mu_oil_fit, w_DME_oil)\n",
    "dμ_oil = w_DME_oil -> polyval(polyder(mu_oil_fit), w_DME_oil)\n",
    "Kval = w_DME_water -> polyval(K_fit, w_DME_water)\n",
    "dKval = w_DME_water -> polyval(polyder(K_fit), w_DME_water)\n",
    "w_oil = w_DME_water -> polyval(Poly([0.0,1.0])*K_fit, w_DME_water)\n",
    "dw_oil = w_DME_water -> polyval(polyder(Poly([0.0,1.0])*K_fit), w_DME_water)\n",
    "KRW  = sw -> krw(sw, krw0, sor, swc, n_w)\n",
    "dKRW = sw -> dkrwdsw(sw, krw0, sor, swc, n_w)\n",
    "KRO  = sw -> kro(sw, kro0, sor, swc, n_o)\n",
    "dKRO = sw -> dkrodsw(sw, kro0, sor, swc, n_o)\n",
    "PC, dPC, d2PC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## physics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: x1 not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: x1 not defined\u001b[39m",
      "",
      "Stacktrace:",
      " [1] \u001b[1mcellEval\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::##93#94, ::JFVM.CellValue{Float64}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m/home/ehsan/.julia/v0.6/JFVM/src/JFVMtools.jl:102\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "# outside the loop: initialization\n",
    "uw       = gradientTerm(p_val) # only for initialization of the water velocity vector\n",
    "k_face   = harmonicMean(k)     # permeability on the cell faces\n",
    "\n",
    "# this whole thing goes inside two loops (Newton and time)\n",
    "\n",
    "\n",
    "\n",
    "∇p0      = gradientTerm(p_val)\n",
    "∇s0      = gradientTerm(sw_val)\n",
    "\n",
    "c_oil         = cellEval(w_oil, c_val) # [mass frac] DME concentration in oil\n",
    "c_oil_face    = arithmeticMean(c_oil)\n",
    "c_face  = arithmeticMean(c_val)\n",
    "sw_face       = upwindMean(sw_val, uw) \n",
    "\n",
    "mu_water_face = faceEval(μ_water, c_face)\n",
    "dμ_water_face = faceEval(dμ_water, c_face)\n",
    "mu_oil_face   = faceEval(μ_oil, c_oil_face)\n",
    "dμ_oil_face   = faceEval(dw_oil, c_face).*faceEval(dμ_oil,c_oil_face)\n",
    "rho_oil_cell  = cellEval(ρ_oil, c_oil)\n",
    "rho_oil_face  = faceEval(ρ_oil, c_oil_face)\n",
    "dρ_oil_cell   = cellEval(dw_oil, c_val).*cellEval(dρ_oil, c_oil)\n",
    "dρ_oil_face   = faceEval(dw_oil, c_face).*faceEval(dρ_oil, c_oil_face)\n",
    "dρc_oil_cell = dρ_oil_cell.*c_oil+cellEval(dw_oil, c_val).*rho_oil_cell\n",
    "dρc_oil_face = dρ_oil_face.*c_oil_face+faceEval(dw_oil, c_face).*rho_oil_face\n",
    "dpc_face      = faceEval(dPC, sw_face)\n",
    "d2pc_face      = faceEval(d2PC, sw_face)\n",
    "\n",
    "dρ_dμ_oil = (dρ_oil_face.*mu_oil_face-dμ_oil_face.*rho_oil_face)./(mu_oil_face.*mu_oil_face)\n",
    "dρc_dμ_oil = (dρc_oil_face.*mu_oil_face-dμ_oil_face.*rho_oil_face.*c_oil_face)./(mu_oil_face.*mu_oil_face)\n",
    "\n",
    "krw_face      = faceEval(KRW,sw_face)\n",
    "kro_face      = faceEval(KRO,sw_face)\n",
    "\n",
    "λ_w_face      = k_face.*krw_face./mu_water_face\n",
    "λ_o_face      = k_face.*kro_face./mu_oil_face\n",
    "\n",
    "dλ_w_face     = k_face.*faceEval(dKRW,sw_face)./mu_water_face\n",
    "dλ_o_face     = k_face.*faceEval(dKRO,sw_face)./mu_oil_face\n",
    "\n",
    "# λ_w_face     = harmonicMean(λ_w) # Wrong!\n",
    "# λ_o_face = harmonicMean(λ_o)     # Wrong!\n",
    "\n",
    "# water mass balance (wmb)\n",
    "M_t_s_wmb, RHS_t_s_wmb = transientTerm(Δsw_init, dt, rho_water*ϕ)\n",
    "M_d_p_wmb = diffusionTerm(rho_water*λ_w_face)\n",
    "M_a_s_wmb = convectionUpwindTerm(rho_water*dλ_w_face.*∇p0)\n",
    "M_a_c_wmb = convectionUpwindTerm(-rho_water*dμ_water_face./mu_water_face.*λ_w_face.*∇p0)\n",
    "\n",
    "# oil mass balance (omb)\n",
    "M_t_s_omb, RHS_t_s_omb = transientTerm(Δsw_init, dt, rho_oil_cell.*ϕ)\n",
    "M_t_c_omb, RHS_t_c_omb = transientTerm(Δc_init, dt, rho_oil_cell.*ϕ)\n",
    "M_d_p_omb = diffusionTerm(rho_oil_face.*λ_o_face)\n",
    "M_a_c_omb = convectionUpwindTerm(dρ_dμ_oil.*k_face.*kro_face.*(∇p0+dpc_face.*∇s0))\n",
    "M_a_s_omb = convectionUpwindTerm(rho_oil_face.*(dλ_o_face.*∇p0+\n",
    "                                (dλ_o_face.*dpc_face+λ_o_face.*d2pc_face).*∇s0))\n",
    "M_d_s_omb = diffusionTerm(rho_oil_face.*λ_o_face.*dpc_face)\n",
    "\n",
    "# DME mass balance (dme)\n",
    "M_t_s_dme, RHS_t_s_dme = transientTerm(Δsw_init, dt, rho_water*c_val-rho_oil_cell.*c_oil)\n",
    "M_t_c_dme, RHS_t_c_dme = transientTerm(Δc_init, dt, rho_water*sw_val+\n",
    "                                       (1.0-sw_val).*dρc_oil_cell)\n",
    "M_d_p_dme  = diffusionTerm(rho_water*c_face.*λ_w_face+rho_oil_face.*c_oil_face.*λ_o_face)\n",
    "M_a_s_dme  = convectionUpwindTerm((rho_water*c_face.*dλ_w_face+\n",
    "                                  rho_oil_face.*c_oil_face.*dλ_o_face).*∇p0+\n",
    "                                  rho_oil_face.*c_oil_face.*\n",
    "                                  (dλ_o_face.*dpc_face+λ_o_face.*d2pc_face).*∇s0)\n",
    "M_a_c_dme  = convectionUpwindTerm((rho_water*(mu_water_face-dμ_water_face.*c_face)\n",
    "                                  ./mu_water_face.*λ_w_face+\n",
    "                                    dρc_dμ_oil.*k_face.*kro_face).*∇p0+\n",
    "                                  dρc_dμ_oil.*k_face.*kro_face.*dpc_face.*∇s0)\n",
    "M_d_s_dme  = diffusionTerm(c_oil_face.*rho_oil_face.*λ_o_face.*dpc_face)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: x1 not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: x1 not defined\u001b[39m",
      "",
      "Stacktrace:",
      " [1] \u001b[1mcellEval\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::##93#94, ::JFVM.CellValue{Float64}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m/home/ehsan/.julia/v0.6/JFVM/src/JFVMtools.jl:102\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "c_oil         = cellEval(w_oil, c_val) # [mass frac] DME concentration in oil"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
